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PREFACE 


A contract  for  an  analytical  study  of  flow  over  toe  cvirves  of  low  spill- 
ways in  gravity  fields  was  negotiated  with  and  awarded  to  the  Curators  of  the 
University  of  Missouri,  Columbia,  Missouri,  by  the  U.  S.  Army  Engineer  Waterways 
Experiment  Station,  The  duration  of  the  contract  was  from  1 September  I966  to 
30  September  I968.  The  Principal  Investigator  was  Dr.  John  J.  Cassidy,  Associ- 
ate Professor  of  Civil  Engineering,  Department  of  Civil  Engineering,  College  of 
Engineering,  University  of  Misso\iri.  Hiis  report  is  the  final,  report  submitted 
under  terms  of  the  contract.  Contract  No.  DACW39-67-C-0001 

The  study  was  financed  by  FI  I967  and  I968  funds  furnished  the  Waterways 
Experiment  Station  by  the  Office,  Chief  of  Engineers,  for  Engineering  Studies 
Program,  ES  8o4,  Development  of  ^^y(iraulic  Design  Criteria, 

Directors  of  the  Waterways  Experiment  Station  during  the  contract  period 
were  COL  John  R.  Oswalt,  Jr.,  CE,  and  COL  Levi  A.  Brown,  CE.  Technical  Director 
was  Mr.  J.  B.  Tiffany.  Mr.  E.  P.  Fortson,  Jr.,  was  Chief,  hydraulics  Division. 
Mr.  F.  B.  Campbell  was  Chief,  hydraulic  Analysis  Branch  during  the  initial 
phase  of  the  contract  study  and  Mr.  E.  B.  Pickett  d\iring  the  final  phase. 

Chief  of  the  Analysis  Section  was  Mr,  R.  G.  Cox. 
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Description  of  The  Problem 

Flow  over  a spillway  is  represented  graphically  in  Figure  1. 

Figure  2 shows  the  portion  of  the  flow  which  was  analyzed  in  this  study. 
The  object  of  the  study  was  to  determine  the  pressure  distribution  along 
the  floor  of  the  toe-curve  and  the  coordinates  of  the  free  surface. 

In  analyzing  the  problem  the  following  assumptions  were  made: 

1.  The  flow  is  inviscid  and  irrotational . 

2.  The  discharge  per  unit  width  Q/L  is  known. 

3.  The  flow  remains  in  contact  with  the  lower  boundary. 

4.  Gravity  acts  in  the  negative  Y direction. 
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The  variables  represented  in  Figures  1 and  2 are  defined  as  follows: 
g ■ Acceleration  due  to  gravity 

H * Total  head 

R » Radius  of  Toe  Curve 

S • Slope  (tan  8) 

Vo  “Approach  velocity 
X,Y  ■ Coordinates 
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*=  Discharge  per  unit  width  of  toe  curve 
PTX,PTY  = X and  Y coordinates  respectively  of  the  point  of  tangency 
PCX.PCY  = X and  Y coordinates  repectively  of  the  point  of  curvature 
XR«YR  = X and  Y coordinates  respectively  of  the  Center  of  curvature 
of  the  toe  curve 
V * Magnitude  of  the  velocity 
3 » Angle  of  inclination  of  downstream  tangent 
a = Angle  of  inclination  of  tangent  at  the  point  of  tangency 
Unless  submerged,  the  spillway  acts  as  a control  section.  In  this 
analysis  it  was  assumed  that  the  spillway  controlled  the  flow  and  that 
the  total  head  H and  the  discharge  per  unit  width  of  toe  curve  q,  could 
be  specified  independently. 

When  formulated  mathematically  the  problem  to  be  solved  is  a 
boundary-value  problem.  The  known  boundary  conditions  are; 

1.  The  pressure  at  the  free  surface  is  zero  (atmospheric). 

2.  At  the  solid  boundary  the  velocity  vector  is  tangent  to  the 
surface. 

In  addition,  boundary  conditions  had  to  be  assumed  at  both  the 
upstream  limit  (Line  A-A')  and  at  the  downstream  limit  (Line  B-B').  At 
the  downstream  end  the  velocity  was  assumed  to  be  constant  (Hydrostatic 
pressure  distribution).  At  the  upstream  end  it  was  assumed  that  the. 
pressure  was  zero  and,  further,  that  a potential  there  was  a straight 
line. 

Mathematical  Statement 

Because  the  pressure  is  zero  along  the  free  surface  the  velocity 
there  is  given  by 

. V - ✓2g(H-Tr  (1) 
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On  the  solid  boundary  the  pressure  head  will.be  given  by 


P 


(2) 


In  order  to  analyze  the  problem  it  is  desirable  to  map  the  physical 
plane  shown  in  Figure  2 into  the  complex-potential  plane  shown  in  Figure  3. 
In  Figure  3 and  4)  represent  the  two-dimensional  stream  function  and 
potential  function  respectively. 


Fig.  3»  Complex-Potential  Plane 

Within  the  complex-potential  plane  the  problem  is  to  determine 
X,  Y and  V as  functions  of  (|>  and  The  velocity  vector  has  a magnitude 
V and  an  inclination  6 which  can  be  expressed  as  the  real  and  imaginary 
parts  of  the  complex  velocity  [1].  The  natural  logarithm  of  the  magnitude 
of  the  velocity  InV  and  the  inclination  e must  satisfy  the  laplace 
equations 


3*6  . 

w 

14  =0 

(3) 

34»^ 

(4) 

within  the  complex-potential  plane  [2].  The  parts  InV  and  e are  also 
related  by  the  Cauchy-Riemann  conditions  [3]. 

(5) 
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In  addition,  connections  between  the  physical-plane  coordinates  and  InV 
and  9 are  provided  by  the  following  integral  equations  valid  only  along 


lines  of  constant  <1)  or  tj/  [2]. 

X 
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X 

II 
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X 
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y2  “ yi  * 

cose 

(9) 

(10) 


The  system  of  equations  (3),  (4),  (5),  (6),  (7),  (8),  (9),  and  (10)  must 
be  solved  in  accordance  with  the  proper  boundary  conditions.  Equation  (1) 
represents  one  of  the  boundary  conditions. 

At  the  upstream  end  the  boundary  condition  must  be  formulated  in 
equation  form.  With  the  pressure  equal  to  zero  along  A-A'  in  Figures 
2 and  3,  the  velocity  must  be  given  by 


V = - = ^2g(K-Tr 


Fig.  4.  Definition  Sketch  for  Upstream  Boundary  Condition 
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Integrating  Eq.  11  with  dr  = dy/cosa 

J = -L-  dy  + C 

3r  cosa 

SgfsiS  [29(H-Y)]3/2  + c (12)- 

where  C is  a constant  of  integration  which  is  evaluated  by  setting 
i|;  = at  Y = PCY. 

At  this  point  it  is  convenient  to  introduce  dimensionless  parameters. 
Equation  (12)  is  made  dimensionless  by  dividing  through  by  . 

Thus,  Eq.  12  becomes 

■ - c,  (13) 

other  variables  expressed  in  dimensionless  form  are 
)L»  B.*  V » q.  _ p , P_ 

Setting  = Cp  at  Y/H  = PCY/H,  Eq.  13  becomes 

= - 2 ([1-?]^''^  + [1-^]^/^)  - Cq  (14) 

3cosa 

where  the  bar  over  the  letter  indicates  the  dimensionless  form  of  the 
variable.  Manipulation  of  Eq.  (14)  yields 

V . [l.?]’/2  . (3cosa(|j,  . ^ (,.j,j^3/2jl/3 

And,  hence 

In?  • I In  (}  - Cp)  + {.\-Kr)^'h 

Equation  (16)  was  used  to  express  In?  In  terms  of  $ along  line  A-A' 

In  Figure  3.  However,  It  was  used  In  the  form 


(15) 

(16) 
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In^  = In  + i In  [{(Vg/V^)^  . 1 } + 1]  (17j 

where  Vg  and  are  the  velocity  at  A and  A'  respectively. 

Numerical  Solution 

Solution  of  the  equations  presented  in  the  previous  section  v/as 
accomplished  in  the  complex-potential  plane  because  the  geometry  is 
regular  there.  An  iterative  procedure  was  used. 

Along  the  solid  boundary  Y = Y(x)  is  known  at  the  outset  but 
Y = Y(4»)  is  required.  A short  FORTRAN  program  was  written  to  generate 
initial  values  for  X(o,il'),  Y((}),ij;),  lnV((|),ii)) , and  This  program 

requires  that  a value  be  assumed  for  the  distance  A-A'  in  Figure  3.  Using 
the  known  values  of  Cg,  H,  and  geometry,  a free  surface  v/as  generated  by 
assuming  that  the  depth  of  flow  was  the  same  at  the  downstream  limit  and 
that  both  and  Y(<p,ip)  varied  linearly  in  the  direction  of  flow. 

This  is  admittedly  a crude  approximation,  but  it  appeared  to  be  satisfactory 
as  a first  approximation.  A copy  of  this  program  is  given  in  the 
appendix.  The  dimensionless  distance  A-A'  is  called  D. 

Using  the  initial  approximations  the  iterative  procedure  was  as 
fol 1 ows : 

1.  Values  of  were  calculated  on  the  free  surface  by 

numerical  differentiation  (dy/dx  = tan'^Gj. 

2.  Values  of  0((J»,i;/)  were  calculated  on  the  solid  boundary  according 
to  the  known  geometry. 

3.  The  0 field  was  settled  using  a numerical  approximation  to  Eq.  3. 

4.  Values  of  InV  were  calculated  on  the  free  surface  using  Eq.  (1). 

5.  Values  of  In?  on  the  solid  boundary  were  calculated  by  numerical 
integration  using  a numerical  approximation  to  Eq.  (6). 
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6.  The  ln7  field  was  settled  using  a numerical  approximation  to 
Eq.  (4). 

7.  New  values  of  X were  computed  on  the  solid  boundary  by  numerically 
integrating  Eq.  (8)  starting  at  the  PC. 

8.  New  values  of  7 on  the  solid  boundary  were  computed  in  accordance 
with  the  X values  obtained  in  step  (7). 

9.  New  values  of  X and  7 on  the  free  surface  were  obtained  by 
numerically  integrating  Eqs.  (7)  and  (9)  respectively,  starting' 
at  the  solid  boundary. 

10.  The  new  values  of  X and  7 were  compared  with  those  used  in  steps 
(1)  and  (2).  If  a significant  change  was  noted  a new  iteration  was 
begun  at  step  (1).  Successive  iterations  were  performed  until 
the  X and  7 values  remained  constant. 

The  iterative  procedure  outlined  in  steps  (1)  through  (10)  was 
programmed  for  processing  by  computer.  The  program  is  included  in  the 
appendix. 

A second  program  was  written  to  calculate  pressures  once  a solution 
was  obtained.  This  program  is  also  contained  in  the  appendix. 

Example  Solution 

A numerical  solution  was  performed  for  the  spillway  toe  curve 
shown  on  plate  34  of  EM  1110-2-1603,  Hydraulic  Design  of  Spillways. 

March  1965.  The  following  table  lists  the  variables  used  in  the  solution, 
their  equivalent  form  in  the  Fortran  language,  and  the  magnitude  assigned 
In  the  example  solution. 


Fortran  Corresponding  Variable  in  Magnitude 

Variable  Figs.  1 and  2 Assigned 
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32.2 

R 

R/H 

2.0 

i 

S 

S 

0.0 

X 

X/H 

Y 

Y/H 

Q 

0.445 

PTX 

PTX/H 

1.893 

PTY 

PTY/H 

-0.676 

PCX 

PCX/H 

0.724 

PCY 

PCY/H 

-0.299 

XR 

XR/H 

1.893 

YR 

YR/H 

1.324 

BETA 

6 

0.0  radians 

ALPH 

a 

-0.6242  radians 

TH 

0 

V 

In  the  computer  programs  TH,  V,  X,  and  Y are  stored  in  arrays  and 
are  dimensioned  as  TH(I.J),  V(I,J),  X(I,J),  and  Y(I,J)  respectively. 

The  subscript  I refers  to  the  row  (1=1  refers  to  the  uppermost  streamline) 
while  J=1  refers  to  the  upstream-most  equipotential  line  — A-A'  in  Fig.  3). 
I and  J have  maximum  values  designated  as  N and  M respectively  for  the 
variables  TH(I,J)  and  V(I,J).  .However,  in  X(I,J)  and  Y(I,J)  I is  either 
1 or  2 while  J has  a maximum  value  equal  to  M.  In  the  example  problem 
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N ■ 9 and  M - 60.  Thus  X(I,J)  and  Y(I,J)  are  2 x 60  arrays  v/hile  V(I,J) 

amd  TH(I,J)  are  both  9 x 60  arrays  in  the  example.  Values  of  TH(I,J)  must 

■» 

be  in  radians  and  TH(I,J)  = 0 corresponds  to  the  positive  X axis. 

Results 

The  pressure>head  distribution  and  free-surface  profile  obtained 
in  the  example  solution  are  shown  in  Fig.  5.  Considerable  difficulty  v/as 
encountered  in  the  vicinity  of  the  upstream  end.  As  can  be  seen  by  the 
computed  free-surface  profile,  some  numerical  error  still  exists  at  the 
upstream  end.  However,  the  pressure  distribution  obtained  is  smooth  and 
was  not  influenced  by  the  irregularity  at  the  free  surface. 

In  the  use  of  the  programmed  solution  it  is  necessary  to  decide 
what  error  will  be  tolerated  in  the  free  surface  coordinates.  In  the 
FORTRAN  program  this  is  labeled  DEL.  In  the  example  solution  DEL  was 
set  equal  to  0.01.  Thus,  the  probable  error  in  X/H  or  Y/H  is  of  the 
order  0.01. 

In  the  relaxation  procedure  used  to  settle  the  TH  and  V fields, 
a tolerable  error  must  also  be  specified.  This  error  has  been  called  ERR. 
Experience  with  the  numerical  procedure  indicates  that  ERR  should  never 
be  greater  than  DEL/10  if  stability  in  the  solution  is  to  be  obtained.  . 

A total  of  33  iterations  were  required  to  obtain  the  solution.  On 
the  IBM  7040  this  required  a total  time  of  6 minutes. 

Numerical  accuracy  is  undoubtedly  closely  related  to  the  curvature 
of  the  streamlines.  Experience  in  this  solution  indicated  that  the 
number  of  streamlines  used  (the  value  of  N)  should  be  such  thst  the 
Increment  in  the  streamfunction  between  streamlines  (Cp/N-1)  should  be 
less  than  0.05.  On  geometry  where  R is  smaller  than  2.0  (as  used  in  the 
example) it  may  be  necessary  to  use  more  streamlines. 
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Use  of  the  Program  Decks 

Schematic  representations  of  the  program  decks  are  contained  in 
the  appendix.  To  use  the  program  for  obtaining  initial  approximations  for 
TH,  V,  X,  and  Y one  must  choose  a value  for  the  dimensionless  distance 
A-A'.  This  is  called  0.  The  program  v/ill  punch  out  cards  containing  the 
required  initial  approximation  to  TH,  V,  Y,  and  X. 

In  the  main  program  a value  for  LIMIT  is  read  in  to  control  the  number 
of  iterations  the  computer  will  perform  {it  may  be  desirable  to  punch 
out  intermediate  values  of  V,  TH,  X,  and  Y because  of  time  limitations). 

A statement  DONE  will  be  typed  out  when  a solution  has  been  obtained. 

The  pressure  computation  deck  uses  the  output  of  the  main  deck  and 
prints  out  pressure  heads  and  free-surface  coordinates. 
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OBJECT  PUOCKAB  IS  HEriG  ft.tElEC  INTO  aTCKAOE. 
SPILLXAV  TCC-CUKVE  PHCSSUi<ES 


PKtSSUiTE  HEAO/H 


FrtCE-SUKFACE  PKCFILE 
X/H  Y/H 


0.72 

-0.00 

0.97 

' 0.0*  1 j ; 

C.77 

0.23 

0.98 

0.0*  ■ : 

C.81 

0.*3 

1.02 

0.03  , ; 

0.87 

0.61 

1.08 

■ O.Ol  ; 

0.92 

0.68 

1.13 

-C.Ol 

0.98 

0.70 

1.18 

-0.0*  . , • 

l.O) 

0.7* 

1.23 

-0.06  ' i ' 

1.09 

0.77 

1.28 

-0.08  ! ■ 

1.18 

0.80 

1.32 

-0.13 

1.21 

0.62 

1.37 

-0.12  1. 

1.27 

0.8* 

l.*2 

-O.l*  ’ 

0.88  T 

l.*7 

,“0.16  ; 

l.AO 

0.86 

1.52 

: -0.17  ' i; 

1.A6 

0.87 

1.57 

-0.19  ‘ ' ’ 

"'n:s2‘ 

0.87 

1.62 

-0.20  ; . 

' 1.88 

0.86 

1.07 

1-0.21  , ii 

‘ 1.6* 

0.88 

1.72 

-0.22 

^r.7o 

0.03 

1.77 

-0.23 

1.77 

0.81 

1.82 

, -0.2* 

> 1.82 

0.78 

1.87 

-0.25  i ' 

1.88 

0.72 

1.92 

-0.25 

1.9* 

0.68 

1.97 

-0.26  i ; 

• 1.99 

' 0.60 

2.02 

, -0.26  1 1 

I ^.08"* 

0.86 

2.07 

; -0.26 

2.10 

0.83 

2.  12 

1 -0.27 

' 2.18 

0.81 

2.17 

-0.27 

— 

0.*9  1 

2.22 

- 2.26 

,o.*7  ; . 

2.27 

.-0.27  j: 

2.31 

c.*t  ■ 

2.32 

-0.28  1.  »• 

■ri2;j6 — 

0.*8 

2.37 

. -0.28 

; 2.*i 

0.** 

2. *2 

: -0.28 

' 2. *6 

'0.*3 

2. *7 

-0.28 

:t^;8i — 

,0.*3 

2.5‘2 

-0.28 

’ : 2.87 

0.*2 

2.57 

’ -0.28  ; ; ; 1 

1 : 2.62 

0.*2 

2.62 

-0.28 

T2T6r 

0.*1 

•2.68 

,:-0.23  i;  il 

, : 2.72 

‘ o.*i  . : 

;2.72 

j : -0.28  , i : 

■ 2.77 

> 0.*l 

■2.70 

11-0.28 

■^“2182 — 
' 2.88 
^2.9J 

o.*o 

o.*o 

o.*o 

2.83 

2.88 

2.93 

' . -0.28  • ; 

; -0.28 

■ -0.28  ■ 

T'2r9a 

1 ; 3.03 
' : 3.08 

o.*o 

'0.*0  ; : 

o.*o  ! • 

“^2.98 

13.03 

'3.09 

1 -0.28  , : ‘ : 

: -0.23  t: 

■ : -0.23  ' ; 

rp3^’i3 — 

: 3.19 
' ; 3.2* 

0.*0  i ; 

, 0.*0  [ ‘ 

: o.*o 

"'3.1* 

i3.19 

13.2* 

■ -0.23 

1 , -0.28  ' : (1 

' ' -o...?o  ( : . ! = 

t ; 3.3* 

0.39 

‘0.39*  ;■  ^ 

3.30 

3.35^ 

; ■ -0.28  ’ „ » : , 

! : -0.28  : ; t 1^ 

Lii  A A « 1 ■ * 

t .......  .!>. .......  .i-c. . . . 1 1 1 n.t.t-1 « ' 1 1 -m  n ftYii « i • » * » * » » --  - 
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program  to  generate  initial  X.Y.V.Tri  FOR  TOE-CURVE  PROGRAM 
DIMENSION  TH)9<60*»V)9>60*tX)2.60»>Y)2>60* 
READ»5tlOO)RtPCXtPCYtPTX.PTY.QtLIMI T 
100  FORMAT (GFlO.St IS) 

REAOtS.lOl )M.NfBETAtStERR.OEL.XR,YH,ALPH 
ICl  FORMAT(2IStAFS*3f3F10.S) 

READtS. 107)0 
107  FORMATIFIO.S) 

MN-M-1 

NM-N-1 

A>NM 

OX-O/A 

X(2.1)*PCX 

Y(2tl )»PCY 

00  1 J*2.M 

1 X)2tJ*«X)2.J-l*,OX 

X) l,l*»X)2f 1*-0*SIN)ALPH* 

X)  1 ,M»-X)2.M* 

A>y.N 

OX«)X) 1 tM*-X) 1 .1«*/A 
DO  2 J-2tMN 

2 X) 1 .J*»X)1 . J-1*,0X 

Y)  ltl*«PCY.O*COS)ALPH* 

00  5 J-2.M 
IF)PTX-X)2.J»)4,3.3 

3 Y)2.J**-SQRT )R*R-)X)2tJ*-XR***2*.YR 
GO  TO  5 

A Y)2,J««PTY.S*)X)2f J*-PTX* 

3 CONTINUE 

Y) ltM*»Y)2.M*.0 
0Y-)Y)1.1»-Y)1.M**/A 
DO  6 J*2tMN 

6 Y)1.J*«Y)1.J-1*-0Y 
00  7 J-ltM 

7 V)1,J*.)AL0G)M-Y)1.J***/2. 

V)Ntl*» )AL0G)H-Y)2»1***/2« 

A«NM 

DV»)V)N.1*-V)1*1*«/A 
00  8 I«2»N 
V) I ,M».V) ItM* 

8 V( I,1).V( I-1.1)^0V 
A«  MN 

OV«(VIN.M)-V(N. 1 ) )/A 
DO  20  J-2.MN 
20  V(NtJ)«V(N. J-1 )+DV 
00  10  J<2«MN 
DV«)V)NtJ*-V) l.J**/A 
DO  10  I>2«NM 

10  V) I ,J*.V) I-1.J*,0V 
00  15  J«1.M 

IF  )X)2tJ*-PTX*12t 12*11 

12  Th)N. J*«aTaN)-)X )2*J*-XR*/ )Y)2*J»-YR** 

GO  TO  13 

11  TH)NfJ*>BETA 

13  00  lA  I«1,NM 

U TH) I . J««TH)N.U* 

15  CONTINUE 

WR1TE)6»103*)X) l.u*»Y) 1»J**X)2*J**Y)2*J**J*1»" 

WRl T£ )6* lOA* ) I V H * J** !•! tN^tJ-l *M* 


WRITE )6i 104O ) ) ThI I tJ** I«l>N*iJ-l .M* 

10%  F0RMAT)9f9.5* 

WRITE  I 7* 102* I ) Thi 2 tJ** I-l»N*tJ-l (M* 
WRlTEl7tl02*> )V) I fJ»,l-l.N*.J-l.M* 

WRITE )7>103»)X) 1 •J»>Y)ltJ*»X)2«J»»Y)2tJ*«J«l»M* 

102  FORMAT)9F8.5* 

103  FORMAT)2F10.3>F13.3iF10*3* 

END 
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C iPILLWAY  TOE-CORVE  ANALYSIS 

01  MENTION  THt9.60).V{9»60).Xi2,60)»Yi2.60) 

C READ  IN  OAT  A 

READIS.lOOlR.PCXtPCYtPTX.PTY.Q.LlMI T 

100  FORMAT (6F lO.S t 15 ) 

REA015>101 )M«NibETAi5.ERRtOEL>XR.YR>ALPH 

101  FORMATI2Ib.AF5.J,3F10.5> 

REAOl  5t  102M  ( V(  I »J)  . 1-1,M  .j.i.M) 

READ! 5. 102) ( ( Th( I « J ) 1 1 - 1 tN ) . J- 1 .M ) 

102  FORMAT (9F8. 5 ) 

READ! 5>1031 (X( 1«J) tYl 1 .J) .X(2«J) •Y(2»J) tJ-ltM) 

103  F0RMAT(2F10.5tF15.5.F10.5) 

c compote  parameters 

MK.«M-2 

ML-M-1 

NM«N-1 

V(N«1  )«Al.Oa<SIKT(  1.0-Y<2>1  ) ) ) 

C NiTIALUE  COUNT  OF  ITERATIONS 
NT  IME»0 

C compute  new  ThETA  on  iOLIO  BOUNDARY 
6 DO  3A  J*liM 

IF(X12. J)-PTX)32#33.33 

32  THlN«J)«ATAN(-(XI2fJ)-XR)/(Y(2»J)-YR) ) 

GO  TO  34 

33  Ih(N.J)«BETA 

34  CONTINUE 

c compute  new  upper  free-surface  thetas 

I«1 

A«NM 

DO  9 J-2.ML 

OYOX»Y(  ItJ-l)*(X(I.J)-X(ItJ't-l))/((X(I.J-l)-X(ItJ)|*(X(IfJ-l)-X(ItJ 
l>l)))-»Y(I>J)*(2.»X(lfJ)-X(I«J-»l)-X(ItJ-l))/(IX(ltJ)-X(liJ't’l))«<X(I 
l.J)-XH  ,J-1  ) ) l-fYU  )*1XI  I tJ)-XU  .0-1  ) )/(  (X(  I .J'fD-XU  .0-1)  )*<X 
111  .0+l)-Xl I .0) ) ) 

9  Th(I.O)»ATAN  (DYOX) 

THI 1 .M) -THl 1 .M-1 ) 

0YDX«Y( l.l)*(2.*X(l.l)-X(1.2)-X(l.3))/nX(l.l)-X(1.2))»IX(l.l)- 
lX(l.3)))^Y(1.2)*(X(l.l)-X<1.3)|/((X<1.2)-Xa.l))*U11.2)- 
1X(1.3)))'.Y(1.3)«(X(1.1)-X(1.2))/((X<1.3>-X(1.1))*IXI1.3)-X(1.2))) 
Th( 1 . 1) «ATAN(0Y0X) 

B» ( TH(N.M)-Th( l.M) )/A 
DO  10  I>1.NM 

10  Th<I.1)«Th<N.1) 

DO  11  1>2.NM 

11  Th(  I .M) -ThI  I-1.M)-fB 
C SETTLE  TH  field 

12  BOX-0. 

DO  14  I-2.NM 
DO  14  0>2.ML 
B^ThII.O) 

Th{I.O)»1Th(1-1.O)-fTH(1.O-1)4TH(I.0'F1)>THU^1.O))/4. 

IF(ABS(8-Th(  1.0)  )-ERR)  l<>>l^tl3 

13  B0X«B0X-»1. 

14  CONTINUE 
IF(BOX) 15.15.12 

C COMPUTE  INCREMENT  IN  STREAM  FUNCTION 

15  A>NM 

00-0/ A 

C COMPUTE  NEW  V ON  FREE  SURFACE 


1 

2 


5 

6 


9 

10 

11 

12 

13 

14 


14B 

149 

151 

152 

153 

154 
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56 

57 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 
124 
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si«o« 

VB-EXP,V(N>1)) 

VT-EXPIV(1«1) ) 

V8VT«(VB/vri**3-l. 

DO  29  J-ltM 

29  Vll*J)-(ALOG(1.0-Ya»J|))/2* 

C compute  new  V at  ENOS 

00  30  I>2iN 
SI-SUO* 

VI  I •!  J-VITtD^ALOGMVBVT^&l/O^l,)/!* 

30  VI  I .M)>Va>M) 

C COMPUTE  NEW  V ON  SOL  10  BOUNOARV 
00  31  J>3>MK 
Jl-J-X 

J3«J'»l 
JA» j»2 
00  31  l«2tN 

I2-I-2 

A-lA.«llHlIl»J3)»THii>U3)-THlll>Jl)-TMl]tUl)| 
b*  Thi  I 1 1 J2  >*Tfil  I > j2i-Tm(  I UJAi-lHi  1 1 JA| 

31  VI 1 1 j )«vi I i>j)-i A«si/Aa. 

VIN«2)«IVINi3l«VlN»ltt/2. 

VI N.MLI ■! VINtMl^VINtMA  » |/2. 

C SETTLE  V riELO 
3 BOX-0, 

00  5 1>2»NH 
DO  5 J«2tML 
b-VI I tJ) 

viitJi-ivii-itJi-fVii^itji^viitj^ij^vii.u-m/A. 
IFIABSIB-VIl tJ) )-ERR)5t5tA 
A BOX>uOX'»l< 

3 CONTINUE 
IF IBOX) 7>7.3 

C COMPUTE  X ON  SOLID  BOUNDARY 
7 A'COSI THiNt 1 I I /EXPI VIN»1  I I 
XBOXbO. 

00  16  J«3«M.2 
C»A 

B»cosi TMiNtu-1 1 1 /E»P'  ;n.j-iii 
A* cos  I Th I N« J| I /EXP I VI N.JI  1 
D«XI 2 .U ) 

bX> I 3. •A* 8. *6-0 ) •00/12* 

XI 2,  j-1  I .<XI 2.  J-1 1*XI2* J-2)*§*»/2. 
BXilA^A.ad+CiaOO/S* 

16  X(2rJ)*IXl2*U-2l«BX*0l/2. 

C compute  new  y ON  solid  BOUNDARY 
00  22  J>2tM 
b«Yl2tU) 

IF IPTX-XI2*J)  ) 19*1S*18 

18  YI2tJ)--SORTlR*R-IXI2*JI-*R»**2>*YR 
GO  TO  20 

19  YI2.J(«PTY^S*IXI2tJ>-PTX) 

20  IF lAbSI YI2*J>-B»-OELl22t22t21 

21  XbOX*XBOX*lt 

22  CONTINUE 

C COMPUTE  Y AND  X ALONG  FREE  SURFACE 
00  26  U>ltM 
P«EXP(V(N»J) I 


F 


I 


i 


0«SIN( Th(N. Jl I/P 
a«cos(Th(n.j) I/p 
OX-0. 

0Y«0. 

Ab«XI Jl 
dA*y ( 1 , J| 

00  23  l«3«Nt2 

C-A 

P»exP<V(A-ltJI  I 

PP«EXP(V(X-2tJ|  I 
B-C0S(  TmU-1»JI  )/P 
A«C0S< IH(X-2»J1 l/PP 
F-0 

E-SlM  iMtX-l  t Jl  I /P 
D*61N( TH<X-2tJ) l/PP 
OX-OX*!  F«-4.*E'*’01  •OQ/3. 

23  OY«OY^lC*A,»b>Al*00/3. 

Y( 1 1 J I- I Yl 1 1 J)*Yl2  •JI>0Y)/2. 

X(  1 .J  I > (X(  1 1 JI-^XI  2 t JI-OXI/2. 

C compare  \e<<  ano  olo  cooroinates 

IF (ABSIXI ltJI-AB)-OEL)AOtAO«Al 
AO  IF (ABS( Y( l.J|-BAI-0Etl26<26*Al 
Al  XBOX*XBOX*i. 

26  CONTINUE 

IF (XbOX127>36t27 

27  NT IME-NT IME*1 
IF(N71ME-LIMIT16.37.37 

36  *RITE(6tlC5l 
105  FORMATiSm  OONEI 

37  I»RIIE(6.107|NT1ME 
107  FORMAT! 151 

•RlTE!7tl02ll(V!ItJ|tI«ltNl,J-l.MI 

»RlTE!7tl02l(nMlI*Jl*l«l»NnJ«l*Ml 

•RlTE!7»103IIXIl*J|.YIl.JltX!2»JliY!2»Jl*J»l*MI 

»RlTEI6»103)<X!l.Jt.Y!ltJ|tXI2»J|»y<2»J)*J»l*M| 

•RlTEIO'lCAlIIVIItJItl-ltNItJ-ltMl 

«RlTEI6*10A)l!TMlI»Jl»I»l»NltJ»l»Ml 

lw4  FORMAT!9F10.5I 
ENO 
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i?rLLi»AY  lOE-CURVE  PRESSURES 

DIMENSION  THl9»60).Vl9.60)tX(2.60»iYl2.60) 
read  I^  data 

nRI tE (61IO6) 

1C6  FORMAT!  29H  SPIlLi^AY  TQE-CURVE  PRESSURES//) 
HEAOIS. lOOlRtPCX.PCYtPlX.PlY.QiLlMl T 

100  FORMAT 16F  10. S.  15  ) 

HEAD! 5. 101 )M,N.BEI A.S.ERK.DEl.XR.YR.ALPH 

101  FORMAI(2I5.AF5.3.3F10.5) 

REA0(5.102) ( (V( I ,J) .l-l.N) .J-l.M) 
REA0(5.102l((TH(Itj|ilal.N)tJ*ltMI 

IC2  FOHMAT(9F8.5) 

READ15.103)(X(l.J).Y(l.J).X(2>J)>Y(2»Ut.U>l.H) 
1C3  FORMAT12F10.5.F15.5.F10.5) 
toSiTE  OUT  headings 
MRIT£<6tlOA) 

luA  FCRMAI(52H  X PRESSURE  hEAD/H  FREE 

WR1TE(6«107) 

.07  FORMAT (52H 

COMPUTE  AND  PRINT  OuT  PRESSURES 
DO  3 J«1.M 

P«1.0-lEXP(VIN.u» )**2)-Yl2.J) 

3 i«R;  7E(6iI0)tX(2«J) 

105  FORMAT(F6.2«F9.2iF20.2tF12.2I 

END 


FREE-SURFACE  PROFILE) 


r/h/zt 


Data  D(2ck 


